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We describe the following new features which signifi- 
cantly enhance the power of the recently developed real-space 
imaginary-time GW scheme (Rieger et al, Comp. Phys. Com- 
mun. 117, 211 (1999)) for the calculation of self-energies and 
related quantities of solids: (i) to fit the smoothly decaying 
time/energy tails of the dynamically screened Coulomb inter- 
action and other quantities to model functions, treating only 
the remaining time/energy region close to zero numerically 
and performing the Fourier transformation from time to en- 
ergy and vice versa by a combination of analytic integration 
of the tails and Gauss-Legendre quadrature of the remaining 
part and (ii) to accelerate the convergence of the band sum 
in the calculation of the Green's function by replacing higher 
unoccupied eigenstates by free electron states (plane waves). 
These improvements make the calculation of larger systems 
(surfaces, clusters, defects etc.) accessible. 

71.15.Th,71.20.-b 



I. INTRODUCTION 

Density-functional calculations provide reliable infor- 
mation about the ground state properties of electron sys- 
tems but give, in principle, no access to the excitation 
spectrum of the system under study. Excitations can 
be described by many-body perturbation theory which 
is, however, at present only computationally feasible for 
real materials in its simplest form, the GW approxima- 
tion of Hedin£rfl The latter gives a comparatively sim- 
ple expression for the self-energy operator, which allows 
the one-particle Green's function of an interacting many- 
electron system to be described in terms of the Green's 
function of a hypothetical non-interacting system with 
an effective potential. The Green's function contains in- 
formation not only about the ground-state density and 
energy but also about the quasiparticle (QP) spectrum. 
The GW approximation has been successfully applied to 
the calculation of QP Jpandstructures of semiconductors 
and other materials for a recent review see Ref. 0. 



The real-space imaginary-time GW method, first pro- 
posed by Rojas et ala and,- in a revised form - described 
in detail by Rieger et alM (we will refer to this paper as 
CPC I in the following) offers a more favourable scaling of 
the computational effort with system-size than conven- 
tional reciprocal-space GW schemesJlJ It substantially 
reduces the computational effort and allows to study 
larger systems than previously possible without resorting 
to further approximations such as plasmon-pole modelscl 
for the energy dependence-of the screened interaction or 
model dielectric functions .til 

The new features outlined in the present paper, partic- 
ularly the new treatment of the (imaginary) time/energy 
dependence, further reduce the computational effort of 
the space-time GW scheme by almost an order of magni- 
tude. This is achieved by fitting the smoothly decaying 
large energy/time tails of all quantities involved in a GW 
calculation to simple model functions and treating the 
remaining time/energy region numerically on a Gauss- 
Legendre grid rather than using an equidistant grid and 
fast Fourier transformations (FFT) from time to energy 
and vice versa. In the new scheme these Fourier transfor- 
mations are performed by a combination of analytic in- 
tegration of the tails and Gauss-Legendre quadrature of 
the remaining part. Another improvement of the method 
concerns the convergence of the calculated Green's func- 
tion with the number of unoccupied eigenstates entering 
the eigenstate (band) sum in the Green's function Eq. 
( |2.ip below. Higher unoccupied eigenstates are approx- 
imated by plane waves. This considerably reduces the 
number of eigenstates and energies which have to be com- 
puted in a density-functional calculation (usually within 
the local density approximation (LDA)) preceding a cal- 
culation of the self-energy with a given accuracy. 

The present paper is organized as follows: first we give 
a brief summary of the real-space imaginary-time GW 
scheme in order to clarify notation in reference to CPC 
I (Section |fi|). Then we describe the new treatment of 
the time/energy dependence (Section III) and the plane- 
wave substitution for accelerating the unoccupied-state 
sum convergence of the Green's function (Section IV). 
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II. SUMMARY OF THE REAL-SPACE 
IMAGINARY-TIME GW METHOD 

In the real-space imaginary- time GW methodi'0 for 
computing electron self-energies and related quantities 
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such as dielectric response functions and quasiparticle 
energies the basic quantities Green's function, dielectric 
response function, dynamically screened Coulomb inter- 
action and self-energy are represented on a real-space grid 
and on the imaginary time axis. In those intermediate 
steps of the calculation where it is computationally more 
efficient to work in reciprocal space and imaginary en- 
ergy we change to the latter representation by means 
of Fourier transforms. The choice of representing the 
time/energy dependence on the imaginary instead of on 
the real axis allows us to deal with smooth, decaying 
quantities which give faster convergence. To obtain the 
self-energy eventually on the real energy axis, we fit a 
model function to the computed self-energy on the imag- 
inary axis, and continue it analytically to the real axis. 
The energy dependence of the dynamically screened in- 
teraction is fully taken into account within the method. 
The computational effort scales quadratically with the 
number of atoms in the unit cell and linearly with the 
number of energjz_points N u used to represent the en- 
ergy dependence.EHl 

First, the zeroth-order Green's function is constructed 
in real space and imaginary time: 
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from the LDA wavefunctions \t n k(r) and eigenvalues e„k- 
Then the RPA irreducible polarizability is formed in real 
space and imaginary time: 

X °(r, r'; tr) = -iG LDA (r, r'; it )Glda(y' , r; -it), (2.2) 

and Fourier transformed to reciprocal space and imagi- 
nary energy and the symmetrised dielectric matrixli3 is 
constructed in reciprocal space, 
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After that the symmetrised dielectric matrix is inverted 
for each k point and each imaginary energy in reciprocal 
space and the screened Coulomb interaction is calculated: 



WdG'(k,Kj) 
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and Fourier transformed to real space and imaginary 
time. From that the self-energy operator 

E(r, r'; it) = iG LD A(r, r'; i T )W{r, v'-it), (2.5) 



and its expectation values (qn|E(ir)|qn) are computed. 
The latter are Fourier transformed to imaginary energy 
and fitted to a model function allowing analytic continua- 
tion onto the real energy axis and evaluation of the quasi- 
particle corrections to the LDA eigenvalues by first-order 
perturbation theory in (S — V^P A ). Since all quantities 
go to zero with increasing r — r' we use a finite cutoff 
region in real space which we call the interaction cell. 
Further details of the method can be found in CPC I. 



III. NEW TREATMENT OF TIME/ENERGY 
DEPENDENCE 

A. motivation and basic idea 

The functions we are dealing with are relatively smooth 
on the imaginary time/energy axis. This allows to em- 
ploy a regular time/energy grid which has the advantage 
that the Fourier transformation from imaginary time to 
imaginary energy and vice versa can be done efficiently 
by fast Fourier transformation (FFTX. However, we still 
need of the order of 100 grid pointstJ for good conver- 
gence (resulting quasiparticle energies converged within 
30 meV with respect to the t/lo grid parameters). This 
point is illustrated by Figure |l| shearing the matrix ele- 
ment of the correlation self-energytil for the uppermost 
valence band of Si at T, calculated (with Au = 0.17 
Hartree) with 30, 60, and 120 FFT grid points (a) on 
the imaginary energy axis and (b) analytically continued 
to real energies. Crucially, the convergence on the imagi- 
nary axis transforms into a convergence of similar quality 
on the real axis upon analytic continuation. Looking at 
the time/energy behaviour of the key quantities, partic- 
ularly those which have to be Fourier transformed such 
as polarizability, screened interaction and the matrix el- 
ements of the self-energy (see Figure ||) we observe that 
they possess nontrivial structure only in the region close 
to it = (iuj = 0) whereas they decay smoothly to zero 
for large imaginary times or energies. The FFT grid has 
to be large enough to take account of the tails (reduce 
aliasing) and at the same time it needs to be sufficiently 
dense to describe the structure in the region close to the 
origin properly. 

This suggests another approach: represent the func- 
tions on a suitably chosen grid in a fixed and compara- 
tively small time/energy interval and fit the large it/ilo 
tails to simple model functions which can be Fourier 
transformed analytically, a method suggested by Blase fl 
al$3 in the context of their earlier mixed-space method.Ej 
For the part handled numerically we choose a Gauss- 
Legendre (GL) grid (linearly transformed from (-1,1) to 
(0,T max ) or (0,uj max ), respectively). This turns out to be 
very efficient since the functions have to be computed for 
a relatively small number of time or energy points only 
and this computation of the functions is much more time- 
consuming than the Fourier transformations themselves 
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which are done by Gaussian quadrature over the numeric, 
cal values and adding the Fourier transform of the tail.tll 
The fit of the tail only needs to be performed whenever 
a quantity has to be Fourier transformed from it to ilo 
or vice versa. The rest of the calculation is restricted 
to the GL grid. Hence the following quantities have 
to be fitted: (1) the polarizability Xgg'O*' * r )> (2) the 
screened Coulomb interaction Wgg' (k, iui) and (3) the 
matrix elements of the correlation part of the self-energy 
(qra|£ c (n-)|qn). 



the polarizability Eq. ( |2.2|) . For that reason we fit the 
large-zr tails of-sach (GG k) element of the polarizabil- 
ity Xgg' (k; ir)cS to a decaying exponential a exp(—br) 
(with b > and for r > only since \° is symmetric 
in t). The two fit parameters a and b are exactly deter- 
mined by fitting two time points: the outermost GL grid 
point and one additional point at 1.3r max . This fitting 
procedure turns out to be very reliable. 

The Fourier transformation from ir to ito is done in 
the following way: 
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FIG. 1. Convergence of the matrix element of the corre- 
lation self-energy (E c (iw)) for the valence-band top V\ hl of 
silicon with respect to uimax with a fixed energy grid spac- 
ing of Au)= 0.17 Har. The top panel shows the calculated 
self-energy on the imaginary axis with the analytically contin- 
ued dependence on real energy shown in the lower panel. The 
lines correspond to u max = 5 Har (dotted), 10 Har (dashed) 
and 20 Har (solid), i. e. 30, 60, and 120 FFT grid points, 
respectively. 

The advantages of this treatment of the time/energy 
dependence are obvious: all quantities have to be 
computed for a much smaller number of imaginary 
times/energies thus saving computational time and re- 
ducing storage requirements while retaining the flexibil- 
ity to accomodate general functional forms of the energy 
dependence and not being restricted to particular forms 
such as plasmon-pole models. 



X°(iUj) = ^2 Pi ix°(in) - a exp(-b\Ti\)]exp(-WjTi) 

oo 

+ dr a exp(—b\T\)exp(—ojjT) 



2 X! Pl [x°(* Tl ) - a exp(-b\n\)] cos(ujjTi) 



2 a 



b 2 + w? 



(3.1) 



with GL grid points tj and ujj, GL weights pi, fit param- 
eters a and 6, and x°(* T ) = ~ *Xgg (k; * T )- 

For a small number of matrix elements Xgg' (k> * T ) 
(typically less than 5% of all matrix elements as long 
as r max is large enough to accomodate all the nontriv- 
ial structure of Xo (*"?")) the large it tails cannot be fit- 
ted to a decaying exponential because they increase or 
change sign. This is only the case for small matrix el- 
ements where the function is already close to zero at 
T m ax anyway. We set x°(l-3T ma;E ) to O.lxo(Tma^) there, 
i. e. choose a reasonable decaying constant which takes 
the (already small) function smoothly to zero. Simply 
setting the matrix element to zero for r > T max would 
render the ensuing fit of VFgg' (k, ioj) unneccessarily dif- 
ficult. Anyhow, T max is a convergence parameter which 
can be varied to check the quality of the results. 



C. Dynamically screened Coulomb interaction 

The large-imaginary-energy tail of the dynamically 
screened interaction Wgg 1 (k, iu) is fitted to the Fourier 
transform of a decaying exponential 



B. Polarizability 

The imaginary-time dependence of the Green's func- 
tion derives from the the decaying exponentials exp(e n ^T ) 
with e„ k < (> 0) for r > (< 0) in Eq. (O). 



The asymptotic behavior at large imaginary times is de- 
termined by the slowest-decaying exponentials and can 
thus be approximated by a single exponential. This 
asymptotic imaginary-time dependence carries over to 



dr a exp(—b\T\)exp(±ic<JT) 



2ab 



b 2 T up 1 /3 2 +lo 2 ' 
(3.2) 



The energy region where W is treated numerically has 
to be large enough to comprise the nontrivial structure 



of W(iuo). We found that 



should be between 3 



and 10 times the plasmon energycj for good convergence. 
We could perform the tail fit along similar lines as that 
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of xoj i- e. subtract the analytic tail function from the 
given imaginary-energy W in (0, u} m ax), Fourier trans- 
form this difference numerically and add the analytically 
given Fourier transform of the function fitted to the tail 
back in. However, for a large number of matrix elements 
W^gg' (k, iuj) the tail fit yields a negative /3 2 because they 
decay more rapidly than 1/uj 2 . In this case the func- 



tion a/ (pr + or) has a pole inside the interval (0, 



which does not allow the analytic Fourier transformation 
to be performed and which makes the numerical Fourier 
transformation of the difference between W and the fit 
function virtually impossible to compute. That is why 
we integrate the analytic tail function from Lo max to oo 



in this case, provided that (3 > 



Part of this 



integral can still be solved analytically whereas the re- 
mainder is treated numerically. The Fourier transform 
W(ir) = — i Wgg' (k, ir) is given by: 

W(lTi) = — Pj W(iliJj) COs(u!jTi) 



3 = 1 



■ du ^2~ l 2 



COs(ijJTi) 



(3.3) 



D. Matrix elements of correlation self-energy 

The matrix elements of the correlation self-energy are 
fitted in a similar way as the polarizability. As the 
asymptotic time dependence of the self-energy is deter- 
mined by that of the Green's function which is 'shorter- 
ranged' (in imaginary time) than W we again fit to a 
decaying exponential a exp(-bT). This time, however, 
we have to perform separate fits on the positive and 
negative half-axis since the self-energy is not symmetric 
in imaginary time. Therefore we obtain two contribu- 
tions to (£ c (zu;)) for positive imaginary energies (with 
E c (ir) = -i(S c (ir))): 

(E c (iuj)) = ^ Pi Pc(«T 4 ) - a + exp(-b+\n\)] 

i=l 

i ■ \ a + 

x exp^—iijJjTi) + ; — 
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Here n and Uj are the GL grid points, pj the GL weights, 
a and /3 2 the fit parameters and ^j max the outermost 
GL grid point. The integral on the right hand side of 
Eq. ( |3.4[ ) is solved numerically using a transformed GL 
grid. It converges rapidly since the integr and is going 
to zero like 1/cj 4 . The second part of Eq. (3.4) is given 
analytically. In this way most Wgg' (k, ico) can be fitted 
except for a small number where /3 2 < —^ 2 nax (this only 
occurs for matrix elements which are small anyway). In 
the latter case we take the correct value at oj max smoothly 
to zero by setting 1 to — 0.9uj 2 lax . Again, the quality of 
the results can be checked by varying io max . 



here n and u)j denote GL grid points, pi GL weigths, 
and a+, a_, b + , and 6_ are the fit parameters. The self- 
energy matrix elements on the negative imaginary-energy 
half-axis are given by (S c (— iui)) = (£ c (iu;))*. 

All matrix elements of S c could be fitted in this way 
for all the systems investigated so far (Si, Ge, GaN). 



E. Tests for bulk silicon and zincblende GaN 

In order to test the stability and accuracy of the tail fit 
and the convergence with the GL grid size we performed 
self-energy calculations for bulk Si (cutoff parameters are 
given in Table | and further details are as in CPC I). 

TABLE I. Cutoff and grid parameters used for the test 
calculations for Si and zincblende GaN, respectively, in section 
[II of the present work. 



parameter 


Si 


GaN 


LDA plane-wave cutoff (in Ry) 


13.5 


50. 


GW plane-wave cutoff a (in Ry) 


19. 


38. 


GW real-space grid 


9x9x9 


12 x 12 x 12 


band cutoff (in Ry) 


10. 


16. 


size of k grid 


4x4x4 


4x4x4 


range T max of time grid b (in a.u.) 


6. 


6. 


size of time (energy) grid b 


15 


15 



a Energy cutoff corresponding to radius of circumscribing 
sphere, see Ref. [lo| 

b This parameter is varied in the tests of Section III . 
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FIG. 2. Matrix elements (S c (ir)) (top panel) and (E c (iw)) 
(lower panel) of the correlation self-energy for the va- 
lence-band top r^/ of silicon. The Gauss-Legendre grid 
points for T max = uj max = 6 a.u. and N T — N u — 15 are 
indicated by filled circles. This figure gives an idea about a 
typical functional form on the imaginary time (energy) axis 
and a typical Gauss-Legendre grid which has to be used to 
Fourier transform this function with good convergence. The 
other quantities which have to be Fourier transformed from 
time to energy or vice versa have similar (or less) structure. 

Figure || gives an idea of the sampling of the matrix 
element of the correlation self-energy for the uppermost 
valence band at T on the imaginary time and imaginary 
energy axis by a GL grid with r max = uj max = 6 a.u. and 
N T = N u = 15. The convergence of two typical quasipar- 
ticle energies with the number of GL grid points and the 
size of the imaginary time/energy region treated explic- 
itly is summarized in Table E3.E3 We chose the lowest con- 
duction state Tf 5 and the valence-band bottom T\. The 
latter is in our experience the slowest-converging quasi- 
particle energy which is most sensitive to details of the 
calculation and cutoff parameters. Obviously, N T = N u 
and T m ax = ^max have to be increased simultaneously 
since the main structure of the functions we are dealing 
with is restricted to a limited region around the origin. 
We observe that 15 GL grid points and T max = uj max = 5 
a.u. are sufficient to converge even the slowest-converging 
valence-band bottom to within 30 meV. The same accu- 
racy is obtained with FFT grids for r max = to m ax = 20 
a.u. and 120 grid points, resulting in QP energies of 3.22 
eV (rf 5 ) and -11.58 eV (T\). 

Our GW calculations for GaN in the zincblcndc struc- 
ture were carried out at the experimental lattice constant 
(a = 8.54 a.u.). The LDA wavefunctions and eigenvalues 
used in the self-energy calculation were obtained from 
a standard plane-wave pseudopotential calculation. The 
Ga 4s and 4p states and the N 2s and 2p states were 
treated as valence states and the soft pseudopotentials of 
Troullier and MartinsEJ were used. The cutoff parame- 



ters are given in Table |. Table III shows the convergence 
of the resulting QP energies at T and X as a function of 
the time/energy GL grid parameters. From these results 
we conclude that the speed of convergence of the QP en- 
ergies of GaN with respect to the imaginary time/energy 
grid is similar to that found for bulk Si. 

TABLE II. Convergence of quasiparticle energies Tf 5 and 
Ti (in eV, top of valence band has been set to zero) for Si 
with respect to Gauss-Legendre grid region r ma x = w m oi (in 
a.u.) and number of grid points 7V T = N u . 











N w 








10 


12 


15 


18 


21 


25 


1 15 


3. 


3.27 


3.28 


3.28 








4. 


3.24 


3.24 


3.24 








5. 




3.23 


3.23 


3.23 






6. 






3.23 


3.23 


3.22 




7. 








3.22 


3.23 


3.22 


rj 


3. 


-11.54 


-11.53 


-11.52 








4. 


-11.58 


-11.52 


-11.52 








5. 




-11.53 


-11.56 


-11.56 






6. 






-11.56 


-11.58 


-11.60 




7. 








-11.63 


-11.58 


-11.59 



TABLE III. Convergence of quasiparticle energies (in eV, 
top of valence band has been set to zero) at the V and X 
point of zincblende GaN with respect to Gauss-Legendre grid 
region T ma x — ^max 

(in a.u.) and number of time/energy grid 

points N T — Nu,. 



N T = N w 


12 


15 


18 


25 


7~max — ^max 


4. 


5. 


6. 


8. 


rf 


-15.88 


-15.95 


-15.93 


-15.94 


n 


3.03 


3.00 


3.00 


2.99 


r 15 


11.56 


11.54 


11.53 


11.52 


xi 


-12.96 


-13.00 


-12.97 


-12.98 


x% 


-6.38 


-6.39 


-6.38 


-6.38 


xi 


-2.66 


-2.66 


-2.66 


-2.66 


xi 


4.43 


4.41 


4.40 


4.39 


X? 


7.91 


7.88 


7.87 


7.86 


xi 


13.23 


13.19 


13.16 


13.14 


xi 


15.32 


15.28 


15.26 


15.24 
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The agreement of our calculated quasiparticle energies 
for Si and GaN with experiment is similar to that of pre- 
vious GW calculations but - in contrast to these earlier 
works - dynamical effects are fully included here. 

Fitting of the tails allows to reduce the time/energy re- 
gion which is treated explicitly by more than a factor of 
three. This saves the same factor in the number of FFT 
grid points. Employing a GL grid enables us to reduce 
the number of points where the functions have to be com- 
puted by another factor between two and three. In total, 
the CPU time, memory and disk space requirements de- 
crease by a factor of seven to eight in comparison with 
the time/energy FFT grid treatment described in CPC 
I. 



IV. PLANE WAVE SUBSTITUTION 

A. motivation and basic idea 

A large number of unoccupied states have to be in- 
cluded in the band sum in the Green's function Eq. (2.1) 
for a proper convergence of the resulting self-energy and 
QP energies. With growing system size it becomes in- 
creasingly difficult to provide such a large number of 
eigenstates by a density-functional calculation since di- 
rect diagonalization may not be computationally feasible 
whereas iterative diagonalization yields only a limited 
number of eigenstates or becomes prohibitively expen- 
sive. Besides that it would be desirable to accelerate the 
convergence of the unoccupied-state sum in order to re- 
duce the computational effort for the calculation of the 
Green's function. 




100 150 
band number 



250 



FIG. 3. Band energies as a function of band number for the 
LDA eigenstates of bulk Si (solid line) and the corresponding 
plane- wave states (dotted line). The latter have been shifted 
upwards by 0.26 Ry (with this shift the energies of LDA and 
PW state number 50 coincide). At higher energies the two 
spectra look remarkably similar (see text). 

On the other hand we expect that the higher the en- 
ergy of an unoccupied state the better it should be ap- 
proximated by a free-electron state (plane wave). This 
is illustrated by Figure || showing the band energy as a 



function of the band number for (LDA) eigenstates of Si 
and the corresponding plane-wave states with wavevec- 
tors K = k + G, G being reciprocal lattice vectors of Si. 
At higher energies the two spectra look remarkably sim- 
ilar if we allow for a constant energy shift between the 
two. Although closer examination shows that that this 
assumption is not fully justified for states with moder- 
ately high energies, which, for symmetry reasons, rather 
resemble linear combinations of several plane waves, the 
sum of all unoccupied states above a certain energy cutoff 
can still be reasonably well described by a corresponding 
sum of plane waves. The aspect of taking proper account 
of the weight of the higher unoccupied states seems to be 
more important than their explicit form. Thus there is 
good reason to expect that the number of unoccupied 
eigenstates which have to be explicitly included into the 
band sum in Eq. (2.1) can be substantially reduced by 
adding a sum of plane waves replacing the omitted higher 
unoccupied states.E3 This is indeed the case as is demon- 
strated by the significant improvement of the convergence 
of the QP energies as a function of the band cutoff upon 
adding a plane-wave (PW) contribution to the Green's 
function, cf. Figures |] and [s] below. 



B. method 



The PW contribution to the Green's function Eq. (2.1) 
takes the following form in real-space: 

AGpw(r,r';jT) = — — exp(iKr)exp(— iKr') 



K 



xezp(--(K a -fcg)) ) (4.1) 

where the K vectors corresponding to energies below the 
lowest PW energy are excluded from the reciprocal-space 
sum. The plane waves are normalized with respect to the 
crystal volume V — VucN^, with being the number 
of k points in the Brillouin zone (BZ) and Vjjc the vol- 
ume of the unit cell. The energies of the PW states are 
measured with respect to an energy zero AE = kg/ 2 
which is determined by adjusting the energy of the high- 
est LDA eigenstate included and the highest PW state 
not inclu ded in the calculation of the Green's function, 
see Eq. below. AGpw can be computed analyt- 

ically by transforming the K sum into an integral by 
Ek v uc/(2tt) 3 J d 3 K, i. e. taking the limit N K -> oo 
and solving the resulting integral. It turns out, how- 
ever, that it is more practical to compute AGpw nu- 
merically instead, even though that is computationally 
slightly more expensive. In this way the contributions 
of the (LDA) eigenstates of the system and the plane 
waves are treated on an equal footing which makes for 
a smoother convergence because of compensation of er- 
rors arising from discretization. Fourier transformation 
of AGpn/(r, r'; it) to reciprocal space results in: 
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AGp W (k,G,G';iT) = 



x ezp(-^(k + G) 2 )<5 G G', (4.2) 



with reciprocal-lattice vectors G,G' of the system with 
|k + G| 2 /2 larger than a given cutoff energy. As Eq. (4.2) 
is diagonal in G it is more efficient to first set up AGpw 
in reciprocal space and then transform it to real space 
before adding it to the Green's function contribution of 
the LDA eigenstates taken into account explicitly. 

In order not to destroy the crystal symmetry of the 
Green's function when adding the PW contribution only 
complete stars of G vectors (groups of plane- wave states 
energetically degenerate at any k in the BZ) must be in- 
cluded into or excluded from AGpw- On the other hand 
the number of G excluded from AGpw should be equal 
to Nb an d s , the number of LDA eigenstates taken into ac- 
count explicitly, which, in turn, has to be determined in 
such a way as to include complete groups of energetically 
degenerate (at any k) LDA eigenstates. These two de- 
mands cannot in general be fulfilled simultaneously. As a 
compromise we calculate AGpw from PW states (k+Gj) 
(ordered with respect to their energy) with weights 



Wi(k) = 



Af 2 (k)-Wi(k) 



i < iVi(k) 
for jVi(k) < i < N 2 (k) 
i > JV 2 (k) 

(4.3) 



where Aq(k) and A 2 (k) are the largest possible/smallest 
possible total number of PW states smaller than/larger 
than Nbands, respectively, both containing complete stars 
of PW states only. This k dependent cutoff ensures 
preservation of both symmetry and number of bands and 
works well in practice. 

In order to account for the difference in the energy 
spectra of LDA eigenstates and plane waves we introduce 
an energy shift 



AE 



\G(Nbands)\ 2 — E cut + E 



VBB 



(4.4) 



between the highest PW state(s) not included and the 
high est LDA eigenstate(s) included in the band sum in 
Eq. (2.1). E cut and Ey b b stand for the energy (at k = 0) 
of the highest LDA state taken into account and the LDA 
valence-band bottom, respectively (both measured with 
respect to the Fermi level which is chosen halfway be- 
tween valence band top and conduction band bottom) 
and \\G (Nb an< is)\ 2 is the energy of plane wave number 
Nbands at k = 0. Taking this energy shift at k = is 
somewhat arbitrary, but it turns out that the resulting 
self-energies and QP energies-are not sensitive to the ex- 
act value of the energy shift. lj 
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FIG. 4. Calculated valence band bottom T\ (top panel, 
the valence band maximum has been set to zero), direct gap 
Ti5 (center), and minimal gap of bulk Si as a function of the 
inverse of Nbands, the number of LDA eigenstates used for 
the calculation of the Green's function. Two sets of data are 
shown, the filled circles (solid lines) refer to calculations with 
LDA eigenstates only whereas the open circles (dashed lines) 
are the results of corresponding calculations where the PW 
contribution was added (see text). The numbers along the 
top axis are the Nbands values used in the respective calcu- 
lations. The number of LDA eigenstates needed to converge 
the quasiparticle energy is significantly reduced by adding the 
PW contribution. 



C. Tests for bulk Si and GaN 

We performed a number of tests in order to assess the 
influence of adding the PW contribution on the conver- 
gence of self-energy and QP energies. The cutoff param- 
eters for our tests for Si are given in Table IV. Further 
calculational details are as in CPC I. Figure^ exhibits 
the convergence of valence band bottom (top panel) , di- 
rect gap at r (center), and minimal gap (bottom) as a 
function of the inverse of Nbands, the number of LDA 
eigenstates used for the calculation of the Green's func- 
tion. Two sets of data are shown in each figure, showing 
results obtained with (open circles) and without (filled 
circles) adding the PW contribution. First of all we ob- 
serve that, as expected, both sets of calculations converge 
to the same answer for \/ Nbands — * 0. However, adding 
the PW contribution dramatically improves convergence; 
the number of eigenstates needed for an accuracy of the 
QP energies of 30 meV is about 70 % smaller when the 
PW contribution is added (see also Table [v| For the 
slowly-converging valence ban d b ottom we find that the 
number of eigenstates in Eq. (|]f]) can be reduced by as 
much as 85 %. 
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FIG. 5. Same as figure ^ for valence band bottom (top 
panel), direct gap F1 (center), and conduction state F1 5 (bot- 
tom) of zincblende GaN. 

Since bulk Si might be perceived as a particularly 
plane-wave like system we also tested the method for 
zincblende GaN. The cutoff parameters are given in Ta- 
ble [V. for further details of the calculation see Section 



[HE above. Figure || showing the band-cutoff depen- 
dence of valence band bottom (top panel), conduction 
band bottom (center), and conduction state T1 5 (bot- 
tom) confirms that the PW contribution to the Green's 
function improves the band-cutoff convergence in the case 
of zincblende GaN, too, although not as much as for bulk 
Si. Adding the PW contribution decreases the number 
of bands needed to converge the quasiparticle energies of 
GaN to an accuracy of 30 meV by 30 to 50 % (see also 
Table [v]). 

TABLE IV. Cutoff and grid parameters used for the test 
calculations for Si and zincblende GaN, respectively, in section 
[ty| of the present work. 



parameter 


Si 


GaN 


LDA plane-wave cutoff (in Ry) 


19. 


50. 


GW plane-wave cutofP (in Ry) 


26. 


50. 


GW real-space grid 


12 x 12 x 12 


15 x 15 x 15 


band cutoff b (in Ry) 


10. 


27. 


size of k grid 


4x4x4 


4x4x4 


range r max of time grid (in a.u.) 


6. 


6. 


size of time (energy) grid 


15 


15 


a Energy cutoff corresponding to 


radius of circumscribing 



sphere, see Ref. [10| 
b This parameter is varied in the tests of Section IV 



TABLE V. Calculated quasiparticle energies at the V and 
X point for Si (in eV) as a function of the number of LDA 
eigenstates Nbands included in the calculation of the Green's 
function. Two sets of results are shown, obtained with in- 
cluding (b) or not including (a) the PW contribution to the 
Green's function (see text). The valence band maximum has 
been set to zero. 



band cutoff [Ry] 

1 ' bands 




6. 
70 


8. 
112 


12. 
186 


14. 

246 


16. 

302 


18. 

500 


Y\ 


(a) 
(b) 


-11.61 
-11.44 


-11.60 
-11.46 


-11.58 
-11.45 


-11.56 
-11.46 


-11.52 
-11.45 


-11.48 


ric/ 
1 15 


(a) 
(b) 


3.20 
3.28 


3.23 
3.27 


3.25 
3.27 


3.26 
3.27 


3.26 
3.27 


3.26 


t 2 


(a) 
(b) 


3.93 
3.92 


3.94 
3.94 


3.95 
3.95 


3.96 
3.96 


3.96 
3.96 


3.96 


XI 


(a) 
(b) 


-7.71 
-7.64 


-7.69 
-7.65 


-7.68 
-7.64 


-7.66 
-7.64 


-7.64 
-7.64 


-7.63 


XI 


(a) 
(b) 


-2.85 
-2.78 


-2.83 
-2.80 


-2.82 
-2.79 


-2.81 
-2.80 


-2.80 
-2.79 


-2.79 


XI 


(a) 
(b) 


1.24 
1.38 


1.27 
1.36 


1.31 
1.36 


1.33 
1.35 


1.35 
1.35 


1.36 


x% 


(a) 
(b) 


10.75 
10.67 


10.75 
10.70 


10.75 
10.70 


10.74 
10.70 


10.72 
10.70 


10.71 



TABLE VI. Same as Table |v| for zincblende GaN. 



band cutoff [Ry] 




12. 


16. 


20. 


28. 


32. 


36. 


N bands 




113 


168 


234 


374 


459 


572 


Y\ 


(a) 


-15.88 


-15.85 


-15.82 


-15.77 


-15.76 


-15.73 




(b) 


-15.56 


-15.64 


-15.66 


-15.68 


-15.69 


-15.70 




(a) 


3.07 


3.12 


3.15 


3.19 


3.20 


3.21 




(b) 


3.38 


3.31 


3.28 


3.24 


3.23 


3.23 


n s 


(a) 


11.58 


11.66 


11.72 


11.79 


11.81 


11.83 




(b) 


12.05 


11.97 


11.92 


11.86 


11.85 


11.84 


xi 


(a) 


-12.93 


-12.91 


-12.90 


-12.87 


-12.86 


-12.85 




(b) 


-12.75 


-12.81 


-12.82 


-12.83 


-12.84 


-12.84 


x v s 


(a) 


-6.39 


-6.34 


-6.30 


-6.26 


-6.24 


-6.22 




(b) 


-6.10 


-6.15 


-6.17 


-6.20 


-6.20 


-6.21 


xt 


(a) 


-2.66 


-2.63 


-2.61 


-2.59 


-2.58 


-2.57 




(b) 


-2.52 


-2.54 


-2.55 


-2.56 


-2.57 


-2.57 


xt 


(a) 


4.46 


4.55 


4.61 


4.67 


4.70 


4.72 




(b) 


4.94 


4.86 


4.81 


4.76 


4.75 


4.74 


X? 


(a) 


7.94 


8.01 


8.06 


8.11 


8.13 


8.15 




(b) 


8.33 


8.26 


8.22 


8.18 


8.17 


8.16 




(a) 


13.27 


13.28 


13.29 


13.30 


13.31 


13.31 




(b) 


13.31 


13.31 


13.31 


13.31 


13.31 


13.31 


xt 


(a) 


15.35 


15.38 


15.40 


15.42 


15.42 


15.43 




(b) 


15.46 


15.45 


15.44 


15.43 


15.43 


15.43 
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It can be concluded from the Fermi energy shifts 
shown in Figure o that adding 



AEp 



E 



QP 



Ep DA 



the PW contribution improves the convergence of abso- 
lute self-energies in qualitatively the same way as that of 
QP energy differences (Figures |^ and ||) . 

In summary we find that the number of LDA eigen- 
states (bands) needed to converge the QP energies within 
30 meV can be considerably reduc ed by including the PW 
contribution described in Section [V B in the calculation 
of the Green's function. 



HI 
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0.05 
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CD 
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0.008 



FIG. 6. Fermi energy shift AE F = E% p - E^ DA for bulk 
Si (top panel) and zincblende GaN (bottom panel) as a func- 
tion of the inverse of Ntands , the number of LDA eigenstates 
used for the calculation of the Green's function, calculated in- 
cluding (open circles, dashed lines) or not including (filled cir- 
cles, solid lines) the plane- wave contribution (see text). The 
plane- wave contribution improves the convergence of absolute 
self-energies in qualitatively the same way as that of the QP 
energy differences shown in figures ^ and [|. 



V. SUMMARY 

In the present work we described two new features 
which significantly enhance the power of the real-space 
imaginary-time GW scheme for the calculation of self- 
energies and related quantities of solids. Fitting the 
smoothly decaying large- imaginary-energy/time tails and 
treating the remaining imaginary energy /time region nu- 
merically on a Gauss-Legendre grid allows to reduce 
the computational time and storage requirements of the 
method by a factor of seven to eight while retaining 
the flexibility to accomodate general functional forms of 
the energy dependence.^ The tail-fitting procedure sug- 
gested in the present work turned out to be accurate 
and reliable. Substituting the contribution of higher un- 
occupied eigenstates to the Green's function Eq. (2.1) 
with a sum of corresponding free-electron states (plane 
waves) acce lerates the convergence of the eigenstate sum 
in Eq. (2.1), thus substantially reducing the number of 



eigenstates and eigenvalues which have to be provided 
by a density-functional calculation preceding the calcula- 
tion of the self-energy and simultaneously decreasing the 
computational effort for the calculation of the Green's 
function itself. 
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